{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Satellite Ground Contacts\n", "\n", "Compute communications contact intervals of \"LandSat 7\" saetllite over a single day" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import satkit as sk\n", "import numpy as np\n", "import plotly.express as px\n", "from datetime import datetime, timedelta\n", "\n", "# The TLE for landsat-7\n", "tle_lines = [\n", " \"0 LANDSAT-7\",\n", " \"1 25682U 99020A 24099.90566066 .00000551 00000-0 12253-3 0 9992\",\n", " \"2 25682 97.8952 129.9471 0001421 108.5441 14.5268 14.60548156329087\"\n", "]\n", "landsat7 = sk.TLE.from_lines(tle_lines)\n", "\n", "# The minimum elevation for a contact\n", "min_elevation_deg = 5\n", "\n", "\n", "# The date to compute ground contacts: April 9, 2024\n", "date = datetime(2024, 4, 9)\n", "# Any array of times representing every second of the day\n", "time_array = np.array([date + timedelta(seconds=x) for x in range(86400)])\n", "\n", "# Get satellite positions in TEME frame (pseudo-inertial) via SGP4\n", "pTEME, _vTEME = sk.sgp4(landsat7, time_array)\n", "\n", "# Get ITRF coordinates (Earth-Fixed) by rotating the position in the TEME frame\n", "# to ITRF frame using the frametransform module\n", "pITRF = np.array([q*x for q,x in zip(sk.frametransform.qteme2itrf(time_array), pTEME)])\n", "\n", "# Setup some ground stations\n", "ground_stations = [\n", " {'name': 'Svalbard', 'lat': 78.2232, 'lon': 15.6267, 'alt': 0},\n", " {'name': 'Alice Springs', 'lat': -23.6980, 'lon': 133.8807, 'alt': 0},\n", " {'name': 'Sioux Falls', 'lat': 43.5446, 'lon': -96.7311, 'alt': 0},\n", "]\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def calc_contacts(ground_station, pITRF, time_array):\n", " \"\"\"\n", " Compute contact times for a single ground station given satellite position in Earth-fixed frame\n", " \"\"\"\n", "\n", " # Create an \"itrfcoord\" object for the ground station\n", " coord = sk.itrfcoord(latitude_deg=ground_station['lat'], longitude_deg=ground_station['lon'], altitude_m=ground_station['alt'])\n", "\n", " # Get the North-East-Down coordinates of the satellite relative to the ground station\n", " # at all times by taking the difference between the satellite position and the ground\n", " # coordinates, then rotating to the \"North-East-Down\" frame relative to the ground station\n", " pNED = np.array([coord.qned2itrf.conj * (x - coord.vector) for x in pITRF])\n", "\n", " # Normalize the NED coordinates\n", " pNED_hat = pNED / np.linalg.norm(pNED, axis=1)[:, None]\n", "\n", " # Find the elevation from the ground station at all times\n", " # This is the arcsine of the \"up\" portion of the NED-hat vector\n", " elevation_deg = np.degrees(np.arcsin(-pNED_hat[:,2]))\n", "\n", " # We can see ground station when elevation is greater than min_elevation_deg\n", " inview_idx = np.argwhere(elevation_deg > min_elevation_deg).flatten().astype(int)\n", "\n", " # split indices into groups of consecutive indices\n", " # This indicates contiguous contacts\n", " inview_idx = np.split(inview_idx, np.where(np.diff(inview_idx) != 1)[0]+1)\n", "\n", " def get_single_contacts(inview_idx):\n", " for cidx in inview_idx:\n", " # cidx are indices to the time array for this contact\n", "\n", " # the North-East-Down position of the satellite relative to\n", " # ground station over the single contact\n", " cpNED = pNED[cidx,:]\n", "\n", " # Compute the range in meters\n", " range = np.linalg.norm(cpNED, axis=1)\n", "\n", " # elevation in degrees over the contact\n", " contact_elevation_deg = elevation_deg[cidx]\n", "\n", " # Heading clockwise from North is arctangent of east/north'\n", " heading_deg = np.degrees(np.arctan2(cpNED[:,1], cpNED[:,0]))\n", "\n", " # Yield a dictionary describing the results\n", " yield {\n", " 'groundstation': ground_station['name'],\n", " 'timearr': time_array[cidx],\n", " 'range_km': range*1.0e-3,\n", " 'elevation_deg': contact_elevation_deg,\n", " 'heading_deg': heading_deg,\n", " 'start': time_array[cidx[0]],\n", " 'end': time_array[cidx[-1]],\n", " 'max_elevation_deg': np.max(contact_elevation_deg),\n", " 'duration': time_array[cidx[-1]] - time_array[cidx[0]]\n", " }\n", " return list(get_single_contacts(inview_idx))\n", "\n", "# Calculate all the contacts\n", "contacts = [calc_contacts(g, pITRF, time_array) for g in ground_stations]\n", "\n", "# Flatten contacts into 1D list\n", "contacts = [item for sublist in contacts for item in sublist]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Convert to pandas dataframe for nice table display\n", "\n", "import pandas as pd\n", "\n", "data = pd.DataFrame(contacts)\n", "data.sort_values(by='start', inplace=True)\n", "data.reset_index(drop=True, inplace=True)\n", "# Get nicer column names for display\n", "data.rename(columns={\"max_elevation_deg\": \"Max Elevation (deg)\",\n", " \"duration\": \"Duration (s)\",\n", " \"start\": \"Start (UTC)\",\n", " \"end\": \"End (UTC)\",\n", " \"groundstation\": \"Ground Station\"}, inplace=True)\n", "data.style \\\n", " .hide(subset=[\"timearr\", \"range_km\", \"elevation_deg\", \"heading_deg\"], axis=1) \\\n", " .format({\"Max Elevation (deg)\": \"{:.1f}\",\n", " \"Start (UTC)\": lambda x: x.strftime('%H:%M:%S'),\n", " \"End (UTC)\": lambda x: x.strftime('%H:%M:%S'),\n", " \"Duration (s)\": lambda x: x.seconds })\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Plot one of the contacts\n", "contact = data.iloc[5]\n", "import plotly.graph_objects as go\n", "from plotly.subplots import make_subplots\n", "\n", "fig = make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.1,\n", " subplot_titles=('Range [km]', 'Elevation [deg]', 'Heading [deg]'))\n", "fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['range_km'], name='Range [km]'), row=1, col=1)\n", "fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['elevation_deg'], name='Elevation [deg]'), row=2, col=1)\n", "fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['heading_deg'], name='Heading [deg]'), row=3, col=1)\n", "fig.update_layout(yaxis = {'title': 'Range (km)'},\n", " yaxis2={'title': 'Elevation [deg]'},\n", " yaxis3={'title': 'Heading [deg]'},\n", " title=f'Landsat 7 to {contact[\"Ground Station\"]} on {contact[\"Start (UTC)\"]}',\n", " width=800,\n", " height=600\n", " )" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.14.0" } }, "nbformat": 4, "nbformat_minor": 2 }